// $Id: dimm.c 484 2014-05-01 00:58:30Z alan.watson.f $

#include "oaxaca-config.h"
#include "oaxaca.h"

#include <time.h>
#include <math.h>

static double dds0[OAXACA_NS_MAX];
static double dds1[OAXACA_NS_MAX];
static double dds2[OAXACA_NS_MAX];

static double interval = 60.0;

static int valid_ns = 10;

static int
valid_dimm(int index)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];
  int ns = detector->ns;
  return (index == OAXACA_DETECTOR_INDEX_WFS) && (ns == valid_ns);
}

static void
dimm_log(void)
{
  FILE *fp = fopen("dimm.csv", "a");
  if (fp == NULL)
    oaxaca_error("dimm: unable to open log file.");

  time_t seconds = time(NULL);
  struct tm tm = *gmtime(&seconds);
  char yyyymmdd[] = "YYYY-MM-DD";
  strftime(yyyymmdd, sizeof(yyyymmdd), "%Y-%m-%d", &tm);
  char hhmmss[] = "hh:mm:ss";
  strftime(hhmmss, sizeof(hhmmss), "%H:%M:%S", &tm);
  fprintf(fp, "\"%s %s\"", yyyymmdd, hhmmss);
  fprintf(fp, ",\"%s\"", yyyymmdd);
  fprintf(fp, ",\"%s\"", hhmmss);

  fprintf(fp, ",%.3f", oaxaca_shm->dimm.interval);
  fprintf(fp, ",%.3f", oaxaca_shm->dimm.r0);
  fprintf(fp, ",%.3f", oaxaca_shm->dimm.fwhm);
  fprintf(fp, ",%.3e", oaxaca_shm->dimm.lambda);
  fprintf(fp, ",%.3e", oaxaca_shm->dimm.K);
  fprintf(fp, ",%.3e", oaxaca_shm->dimm.sigma);
  fprintf(fp, "\n");
  
  fclose(fp);  
}

void
oaxaca_dimm_start(int index)
{
  if (!valid_dimm(index))
    return;
    
  oaxaca_shm->dimm.interval = interval;
  
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];
  int ns = detector->ns;

  for (int i = 0; i < ns; ++i)
    dds0[i] = dds1[i] = dds2[i] = 0.0;

  oaxaca_shm->dimm.sigma  = 0.0;
  oaxaca_shm->dimm.K      = 0.0;
  oaxaca_shm->dimm.r0     = 0.0;
  oaxaca_shm->dimm.lambda = 0.0;
  oaxaca_shm->dimm.fwhm   = 0.0;
}

void
oaxaca_dimm_next(int index)
{
  if (!valid_dimm(index))
    return;

  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];
  int ns = detector->ns;
  int na = OAXACA_NA_MAX;

  double ds[OAXACA_NS_MAX];

  for (int i = 0; i < ns; ++i) {
    ds[i] = detector->s[i] - detector->r[i];
    for (int j = 0; j < na; ++j) {
      ds[i] += detector->i_matrix[j][i] * oaxaca_shm->wfc.a[j];
    }
  }
  

#if 0

  // 5x5 ring

  // We measure longitudinal displacements between pairs of spots in the
  // same row or column with b = 4. The spots must be defined in the
  // correct order, from left to right first, and from top to bottom
  // second. We have six pairs of spots.

  double b = 4.0;
  int np = 6;
  double dds[np];
  dds[0] = ds[ 9 * 2 + 0] - ds[ 0 * 2 + 0];
  dds[1] = ds[10 * 2 + 0] - ds[ 1 * 2 + 0];
  dds[2] = ds[11 * 2 + 0] - ds[ 2 * 2 + 0];
  dds[3] = ds[ 4 * 2 + 1] - ds[ 3 * 2 + 1];
  dds[4] = ds[ 6 * 2 + 1] - ds[ 5 * 2 + 1];
  dds[5] = ds[ 8 * 2 + 1] - ds[ 7 * 2 + 1];

#elif 0

  // 2x2 square

  // We measure longitudinal displacements between pairs of spots in the
  // same row or column with b = 1. The spots must be defined in the
  // correct order, from left to right first, and from top to bottom
  // second. We have four pairs of spots.

  double b = 1.0;
  int np = 4;
  double dds[np];
  dds[0] = ds[2 * 2 + 0] - ds[ 0 * 2 + 0];
  dds[1] = ds[3 * 2 + 0] - ds[ 1 * 2 + 0];
  dds[2] = ds[2 * 2 + 1] - ds[ 0 * 2 + 1];
  dds[3] = ds[3 * 2 + 1] - ds[ 1 * 2 + 1];

#else

  // 3x3 cross

  // We measure longitudinal displacements between pairs of spots in the
  // same row or column with b = 2. The spots must be defined in the
  // correct order, from left to right first, and from top to bottom
  // second. We have two pairs of spots.

  double b = 2.0;
  int np = 2;
  double dds[np];
  dds[0] = ds[4 * 2 + 0] - ds[ 0 * 2 + 0];
  dds[1] = ds[3 * 2 + 0] - ds[ 1 * 2 + 0];

#endif


  for (int i = 0; i < np; ++i) {
    dds0[i] += 1.0;
    dds1[i] += dds[i];
    dds2[i] += dds[i] * dds[i];
  }

  if (detector->frame > 0 && detector->frame % (unsigned long) (interval * detector->framefrequency) == 0) {
  
    double sumsq = 0.0;
    for (int i = 0; i < np; ++i) {
      double variance = (dds0[i] * dds2[i] - dds1[i] * dds1[i]) / dds0[i] / (dds0[i] - 1.0);
      sumsq += variance;
    }
    double sigma = sqrt(sumsq / np);
    oaxaca_debug("dimm: sigma = %.3f pixels (raw)", sigma);
    
    oaxaca_debug("dimm: ESHM is %.2e.", detector->eshm);

    if (detector->eshm > 0) {

      double sigmanoise = sqrt(2) * 57.294 * pow(detector->eshm, -0.8748);
      oaxaca_debug("dimm: sigma = %.3f pixels (noise)", sigmanoise);
    
      double sigmacorrected = (sigmanoise > sigma) ? 0.0 : sqrt(sigma * sigma - sigmanoise * sigmanoise);
      oaxaca_debug("dimm: sigma = %.3f pixels (corrected)", sigmacorrected);
  
      double D = 0.13;
      double lambda = 750e-9;
      double K = 0.364 * (1.0 - 0.532 * pow(b, -1.0 / 3.0) - 0.024 * pow(b, - 7.0 / 3.0));
      oaxaca_debug("dimm: K is %.3f.", K);

      double pi = 4.0 * atan(1.0);
      double arcsec = pi / (180.0 * 60.0 * 60.0);

      double r0 = pow(K, 0.6) * pow(lambda / (sigmacorrected * detector->scale * arcsec), 1.2) / pow (D, 0.2);
      oaxaca_info("dimm: r0 is %.1f cm at %.0f nm.", r0 * 100, lambda * 1e9);

      double fwhm = (r0 == 0.0) ? 0.0 : 0.98 * lambda / r0 / arcsec;
      oaxaca_info("dimm: seeing FWHM is %.2f arcsec.", fwhm);
    
      oaxaca_shm->dimm.sigma  = sigmacorrected * detector->scale;
      oaxaca_shm->dimm.K      = K;
      oaxaca_shm->dimm.r0     = r0;
      oaxaca_shm->dimm.lambda = lambda;
      oaxaca_shm->dimm.fwhm   = fwhm;

      dimm_log();

      for (int i = 0; i < ns; ++i)
        dds0[i] = dds1[i] = dds2[i] = 0.0;

    }

  }

}
